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We present a detailed analysis of the modulational instability of the zone-boundary mode for one 
and higher-dimensional Fermi-Pasta-Ulam (FPU) lattices. Following this instability, a process of 
relaxation to equipartition takes place, which we have called the Anti-FPU problem because the 
energy is initially fed into the highest frequency part of the spectrum, at variance with the original 
FPU problem (low frequency excitations of the lattice). This process leads to the formation of chaotic 
breathers in both one and two dimensions. Finally, the system relaxes to energy equipartition on 
time scales which increase as the energy density is decreased. We show that breathers formed 
when cooling the lattice at the edges, starting from a random initial state, bear strong qualitative 
similarities with chaotic breathers. 

Keywords:Modulational instability, Localized modes, energy localization, energy equipartition, 
chaotic dynamics. 

PACS numbers: 05. 45. -a Nonlinear dynamics and nonlinear dynamical systems 95.10.Fh Chaotic dynamics 
63.20.Pw Localized modes 

Several nonlinear physical systems exhibit modulational instability, which is a self-induced mod- 
ulation of the steady state resulting from a balance between nonlinear and dispersive effects. This 
phenomenon has been studied in a large variety of physical contexts: fluid dynamics, nonlinear optics 
and plasma physics. The Fermi-Pasta-Ulam (FPU) lattice is an extremely well suited model system 
to study this process. Both the triggering of the instability and its further evolution can be studied in 
detail, exciting initially high-frequency modes. The original FPU problem was casted instead in the 
context of long wavelengths. This is why we call the process we analyze in this paper, the Anti FPU 
problem because of the analogy with the seminal FPU numerical simulation. At variance with the 
appearance of (m)KdV-solitons in the FPU original problem, in this process the pathway to equipar- 
tition leads to the creation of localized objects that are chaotic breathers. Similar localized structures 
emerge when cooling the lattice at the edges, starting from thermalized initial states. 

I. INTRODUCTION 

In 1955, reporting about one of the first numerical simulations, Fermi, Pasta and Ulam (FPU) remarked that it 
was . . . very hard to observe the rate of "thermalization" or mixing ... in a nonlinear one-dimensional lattice in which 
the energy was initially fed into the lowest frequency mode. Even if the understanding of this problem advanced 
significantly afterwards 0, 0, several issues are still far from being clarified. In most cases, the evolution towards 
energy equipartition among linear modes has been checked considering an initial condition where all the energy of the 
system is concentrated in a small packet of modes centered around some low frequency. 

Beginning with the pioneering paper of Zabusky and Deem Q, the opposite case in which the energy is put into 
a high frequency mode has been also analyzed. In this early paper, the zone-boundary mode was excited with an 
added spatial modulation for the onc-dimcnsional cv-FPU model (quadratic nonlinearity in the equations of motion). 
Here, we will study the time-evolution of this mode without any spatial modulation for the /3-FPU model (cubic 
nonlinearity in the equations of motion) and some higher-order nonlinearities. Moreover, we will extend the study to 
higher dimensional lattices. Since the energy is fed into the opposite side of the linear spectrum, we call this problem 
the Anti-FPU problem. 

In a paper by Bundinsky and Bountis [5j, the zone-boundary mode solution of the one-dimensional FPU lattice 
was found to be unstable above an energy threshold E c which scales like 1/N, where N is the number of oscillators. 
This result was later and independently confirmed by Flach [|| and Poggi et al. @|, who also obtained the correct 
factor in the large TV-limit. These results were obtained by a direct linear stability analysis around the periodic orbit 
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corresponding to the zone-boundary mode. Similar methods have been recently applied to other modes and other 
FPU-like potentials by Chechin et al HE! and Rink [ljj. 

A formula for E c , valid for all N, has been obtained in Refs. [TH. IT2I Kh\ IT3| in the rotating wave approximation, 
and will be also discussed in this paper. Associated with this instability is the calculation of the growth rates of mode 
amplitudes. The appropriate approach for Klein-Gordon lattices was first introduced by Kivshar and Peyrard |l5j| . 
following an analogy with the Benjamin- Feir instability in fluid mechanics 

Previously, a completly different approach to describe this instability was introduced by Zakharov and Shabat [l 7| , 
studying the associated Nonlinear Schrodinger equation in the continuum limit. A value for the energy threshold was 
obtained in Ref.01 m the continuum limit. The full derivation starting from the FPU equation of motions was then 
independently obtained by Berman and Kolovskii [l9| in the so-called "narrow-packet" approximation. 

Only very recently the study of what happens after modulational instability develops has been performed for 
Klein-Gordon [23 and FPU-lattices[H |2l|. From these analyses it turned out that these high-frequency initial 
conditions lead to a completely new dynamical behavior in the transient time preceeding the final ener gy e quipartition. 
In particular, the main discovery has been the presence on the lattice of sharp localized modes j2{J, l21| . These 
latter papers were the first to make the connection between energy relaxation and intrinsic localized modes or 
breathers I23l 24 l25l |26( . Later on, a careful numerical and theoretical study of the dynamics of a /3-FPU model was 
performed |27j . It has been shown that moving breathers play a relevant role in the transient dynamics and that, 
contrary to exact breathers, which are periodic solutions, these have a chaotic evolution. This is why they have been 
called chaotic breathers. Following these studies, Lepri and Kosevich and Lichtenberg and coworkers p^.l30| have 
further characterized the scaling laws of relaxation times using also continuum limit equations. 

On the other hand, studies of the asymptotic state of the FPU lattice d yna mics when energy is extracted from the 
boundaries have revealed the persistence of localized modes |U H3, 13113 • Already some of these authors [U 
have discussed the similarities of these modes with chaotic breathers. In this paper, we will further study this 
connection. 

Most of the previous studies are for one-dimensional lattices. Here, we will derive modulational instability thresholds 
also for higher dimensional lattices and we will report on a study of chaotic breathers formation in two-dimensional 
FPU lattices. 

We have organized the paper in the following way. In Section [HI the modulational instability of zone-boundary 
modes on the lattice is discussed, beginning with the one-dimensional case, followed by the two-dimensional and 
higher dimensional cases and finishing with the continuum Nonlinear Schrodinger approach. Section IIIII deals with 
the mechanisms of creation of chaotic breathers in one and two dimensions. Finally, in Section IIVI we discuss the 
relation with numerical experiments performed when the lattice is cooled at the edges. Some final remarks and 
conclusions are reported in Section Ivl 



II. MODULATIONAL INSTABILITY 



A. The one-dimensional case 



We will discuss in this section modulational instability for the one-dimensional FPU lattice, where the linear coupling 
is corrected by a (2p + l)-th order nonlinearity, with p a positive integer. Denoting by u n (t) the relative displacement 
of the n-th particle from its equilibrium position, the equations of motion are 

u n = u n+ i + it n _i - 2u n + (u n +i - u n ) 2p+1 - (u n - u n -i) 2p+1 . (1) 

We adopt a lattice of N particles and we choose periodic boundary conditions. For the sake of simplicity, we first 
report on the analysis for p — 1 (i.e. the /3-FPU model) and then we generalize the results to any p-value. 

Due to periodic boundary conditions, the normal modes associated to the linear part of Eq. Q are plane waves of 
the form 

u n (t) = |(e M »W+e- lfl -W) (2) 

where n (t) = qn — ujt and q = 2-kU/N (k = —N/2, . . . ,N/2). The dispersion relation of nonlinear phonons in 
the rotating wave approximation ^| is a - ,2 ( ( ?) = 4(1 + a) sin 2 (q / 2) , where a = 3a 2 sin 2 (q/2) takes into account 
the nonlinearity. Modulational instability of such a plane wave is investigated by studying the linearized equation 
associated with the envelope of the carrier wave J2J). Therefore, one introduces infinitesimal perturbations in the 
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amplitude and phase and looks for solutions of the form 



i{t) = l[l + bn (f)] e *[«»W+^(*)] + £[l +6n ( t) ] e -i[«n(t)-Hfr„(t)] 

= a[l + b n (t)] cos[qn — ut + ip n (t)], 



(3) 



where 6„ and ^ n are reals and assumed to be small in comparison with the parameters of the carrier wave. Substituting 
Eq. iJSJ into the equations of motion and keeping the second derivative, we obtain for the real and imaginary part of 
the secular term e l ^ qn ~ ut ' the following equations 



u> b r , 



(1 + 2a) [cosq (b n+1 + b n -i) - 2b n ] 

a (b n+1 + b n -i ~ 2b n cosq) - (1 + 2a) smq(ip n +i - 4>n-i) 



-u) 2 i/j n - 2ujb n + ip n = (1 + 2a) [cosq (ip n +i + tpn-i) - 2ip n ] 



+ (1 + 2a)sing(6 n+ i - b n -i) + a (ip„+i + t/j n -i - 2^ n cosq) 



(4) 



(5) 



Further assuming b n — bo e 1 ^ 71 + c.c. and tp n — ipo e l ^ n at ' + c.c. we obtain the two following equations for 
the secular term e l (Q n - nt ) 



b 

ipo 



tt 2 + u 2 + 2(1 + 2a) (cos q cos Q - 1) - 2a(cosQ - cosq) - 2itp [ojft + (1 + 2a) sin q sin Q] = (6) 
if + uj 2 + 2(1 + 2a) (cos q cos Q - 1) + 2a(cosQ - cosq) + 2ib [uSl + (1 + 2a) sing sin Q] = 0. (7) 



In the case of Klein-Gordon type equations 0,H(j, one neglects the second order derivatives in Eqs. Q-©- This 
can be justified by the existence of a gap in the dispersion relation for q = 0, which allows to neglect fl 2 with respect 
to lo 2 . In the FPU case, this approximation is worse, especially for long wavelengths, because there is no gap. 

Non trivial solutions of Eqs. ©-0 can be found only if the Cramer's determinant vanishes, i.e. if the following 
equation is fulfilled: 



(CI + u) 2 - 4(1 + 2a) sin 2 



q + Q 



(Q - uj) 2 - 4(1 + 2a) sin 2 
4a 2 (cos Q — cos q) 2 . 



q-Q 



(8) 



This equation admits four different solutions when the wavevectors q of the unperturbed wave and Q of the 
perturbation are fixed. If one of the solutions is complex, an instability of one of the modes (g± Q) is present, with a 
growth rate given by the imaginary part of the solution. Using this method, one can derive the instability threshold 
amplitude for any wavenumber. A trivial example is the q = case, for which we obtain Q = ±sin(<5/2), which 
proves that the zero mode solution is stable. This mode is present due to the invariance of the equations of motion © 
with respect to the translation u n — > u n + const and, as expected, is completely decoupled from the others. 

A first interesting case is q = it. One can easily see that Eq. (JSJ admits two real and two complex conjugate 
imaginary solutions if and only if 



9 Q 1 + a 

cos^ — > . 

2 1 + 3a 



(9) 



This formula was first obtained by Sandusky and Page (Eq. (22) in Ref. [13| ) using the rotating wave approximation. 
The first mode to become unstable when increasing the amplitude a corresponds to the wavenumber Q — 2w/N. 
Therefore, the critical amplitude a c above which the q = 7r-mode looses stability is 



sin 2 (tt/JV) 



, 3 [3 cos 2 (tt/JV) - 1] 
This formula is valid for all even values of N and its large JV-limit is 



1/2 



%/6JV 



O 



1 

JV3 



(10) 



(11) 



In Fig.^ we show its extremely good agreement with the critical amplitude determined from numerical simulations. 
It is interesting to emphasize that the analytical formula (jl()|) diverges for N — 2, predicting that the 7r-mode is 
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stable for all amplitudes in this smallest lattice. This is in agreement with the Mathieu equation analysis (see Ref. 
p. 265). 
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FIG. 1: Modulational instability threshold amplitude for the 7r-mode versus the number of particles in the one-dimensional 
FPU lattice. The solid line corresponds to the analytical formula (I10L the dashed line to its large JV-estimate (I11H and the 
diamonds are obtained from numerical simulations. 



It is also interesting, to express this result in terms of the total energy to compare with what has been obtained 
using other methods |5l Irl ItI H7I Il9| . Since for the 7r-mode the energy is given by E — N(2a 2 + 4a 4 ), we obtain the 
critical energy 



E r = 



2N 



tt\ 7 cos 2 (n/N) - 1 
NJ [3 cos 2 (n/N)-lf 



For large N, we get 



TT 2 f 1 



(12) 



(13) 



This asymptotic behavior is the same as the one obtained using the narrow packet approximation in the context of the 
Nonlinear Schrodinger equation by Berman and Kolovskii (Eq. (4.1) in Ref. 19]). The correct scaling behavior with 
N of the critical energy has been also obtained by Bundinsky and Bountis (Eq. (2.22) in Ref. 0) by a direct linear 
stability analysis of the 7r-mode. The correct formula, using this latter method, has been independently obtained by 
Flach (Eq. (3.20)) in Ref. @) and Poggi and Ruffo (p. 267 of Ref. 0). Recently, the A^-scaling of formula lO has 
been confirmed using a different numerical method and, interestingly, it holds also for the 27r/3 and 7r/2 modes 36] . 

This critical energy is also very close to the Chirikov "stochasticity threshold" energy obtained by the resonance 
overlap criterion for the zone boundary mode[37|. The stochasticity threshold phenomenon has been thoroughly 
studied for long wavelength initial conditions, and it has been clarified that it corresponds to a change in the scaling 
law of the largest Lyapunov exponent 38]. We will show in Section ITTT1 that above the modulational instability critical 
energy for the 7r-mode one reaches asymptotically a chaotic state with a positive Lyapunov exponent, consistently 
with Chirikov's result. 

The above results can be generalized to nonlinearities of 2p + 1 order in the equations of motion . We limit the 
analysis to the 7r-mode, for which the instability condition takes the form 



2 Q 1 + a 
cos — > — -— , 

2 l + (2p+l)a' 



(14) 
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where 



(2P+1)! 2p 

a = — — — a . 

p!(p + l)! 

Hence the critical amplitude above which the 7r-mode is unstable is 

2/ „n nV2 



p!(p+l)! sin 2 (tt/TV) 



(2p + 1)! [(2p + 1) cos 2 (tt/TV) - 1] 



leading to the large TV scaling 



a c 



jV-i/p 



(15) 



(16) 



(17) 
(18) 



This scaling also corresponds to the one found in Ref. |39| when discussing tangent bifurcations of band edge plane 
waves in relation with energy thresholds for discrete breathers. Their "detuning exponent" z has a direct connection 
with the nonlinearity exponent p — z/2. We will see in Section lll Bl that this analogy extends also to higher dimensions. 

For fixed TV, a c is an increasing function of the power of the coupling potential with the asymptotic limit linip^oo a c = 
0.5. Therefore, in the hard potential limit the critical energy for the 7r-mode increases proportionally to TV. The fact 
that we find a higher energy region where the system is chaotic is not in contradiction with the integrability of the 
one-dimensional system of hard rods 40] , because in the present case we have also a harmonic contribution at small 
distances. 

For the FPU-a model (quadratic nonlinearity in the equations of motion) , the 7r-mode is also an exact solution which 
becomes unstable at some critical amplitude which, contrary to the case of the FPU-/3 model, is TV-independent |sl : 
which means that the critical energy is proportional to TV and then that 7r-mode can be stable in some low energy 
density limit also in the thermodynamic limit. 

It has also been realized 0, II Efl El El E3 that group of modes form sets which are invariant under the dynamics. 
The stability analysis 0, of pair of modes has shown a complex dependence on their relative amplitudes. The 
existence of such invariant manifolds has also allowed to construct Birkhoff-Gustavson normal forms for the FPU 
model, paving the way to KAM theory EH . 



B. Higher dimensions 



In this section, we will first discuss modulational instability of the two-dimensional FPU model. The method 
presented in section lll Al can be easily extended and the global physical scenario is preserved. However, the scaling 
with TV of the critical amplitude changes in such a way to make critical energy constant, in agreement with the analysis 
of Ref. H3. 

The masses lie on a two-dimensional square lattice with unitary spacing in the (x, y) plane. We consider a small 
relative displacement u„, m (n, to € [1, TV]) in the vertical z-direction. Already with a harmonic potential, if the spring 
length at equilibrium is not unitary, the series expansion in u n%m of the potential contains all even powers. We retain 
only the first two terms of this series expansion. After an appropriate rescaling of time and displacements to eliminate 
mass and spring constant values, one gets the following adimensionallizcd equations of motions 



^n,m ^n+l,m ~t~ l^n— l,m ^n,m+l ^n,m — 1 4zZ n?n 

~t~ (^n+l,m ^n,m) ~l~ {^n— l,m ^n,m) (^n,m+l ^n,m) (^n.m — 1 ^n,m) 

Considering periodic boundary conditions, plane waves solutions have the form 

u n ,m = a cos (q x n + q y m — ojt) . 
In the rotating wave approximation^^, one immediately obtains the dispersion relation 



lu =4 sin — +4 sin -£ + i2a 



. 4 Qx 
sin — 
2 



• 4 1y 

sin — 
2 



(19) 
(20) 

(21) 
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which becomes exact for the zone-boundary mode (q x , q y ) = (it, tt), 

< w = 8(l + 3o a ). (22) 

In order to study the stability of the zone-boundary mode, we adopt a slightly different approach. Namely, we consider 
the perturbed relative displacement field of the form 

u n , m = (| + 6„, m ) e *(™+™-"'.-*) + c.c, (23) 

where 6„ jm is complex. This approach turns out to be equivalent to the one of Section Til Al in the linear limit. 
Substituting this perturbed displacement field in Eqs. I|19() . we obtain 

[1 + 2a] [b n+ i^ m + b n -i l7n + frn/m+l + &n,?n-l + ^>n,m] 

-a [bn+l,m + K-l,m + K,m+1 + K,m-1 + 4 K,m] = ~K, m + 2iw 7r , 7r & n ,m + c4,-A,m, (24) 

where a = 3a 2 . Looking for solutions of the form 

we arrive at the following set of linear algebraic equations for the complex constants A and B 

[(ft + w^) 2 -8(1 + 2a) A] A + 8aAB = (26) 
8aAA + [(fl - oj-K.-n) 2 - 8(1 + 2a) A] B = 0, (27) 

where 2A = cos 2 (Qa;/2) + cos 2 (Q a /2). As for the one-dimensional case, we require that the determinant of this linear 
system in A and B vanishes, which leads to the following condition 

[(CI + lj^^) 2 - 8A(1 + 2a)] [(CI - Uit^f - 8A(1 + 2a)] = 64a 2 A 2 . (28) 

This equation admits two real and two complex conjugated imaginary solutions in CI if 

which is the analogous for two dimensions of condition @. One can achieve the minimal nonzero value of the r.h.s. 
of the above expression choosing Q x = 0, Q y = 2tt/N, which leads to the following result for the critical amplitude 



sin 2 (7r/A0 > 1 " 



3 [3 cos 2 (tt/A0 + 1] 



(30) 



12 AT V^ 3 



Its large N limit is 



This prediction is compared with numerical data in Fig. [3 The agreement is good for all values of N. 

Since the relation between energy and amplitude is now E = 2N 2 (2a 2 + 4a 4 ), we obtain the critical energy in the 
large AT-limit as 

Ee = 4 + o(±)- (32) 



3 V^ 2 , 

This shows that the critical energy is now constant in the thermodynamic limit, which agrees with the remark of 
Ref. [3^ about the existence of a minimal energy for breathers formation [2(|. 

The results of this Section can be easily extended to any dimension d. Repeating the same argument, we arrive at 
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FIG. 2: Modulation instability threshold for the (tt,-k) mode versus number of oscillators in two dimensional array. The solid 
line is given by the corresponds to the estimate obtained from the Nonlinear Schrodinger equation in large TV limit (see formula 
and the diamonds are results of numerical simulations. 

the following estimates for the critical amplitude and energy in the large TV limit 

71 1 o(±) (33) 



TV 2 

E c = —N d - 2 + O (7V rf - 4 ) . (34) 
o 

This means that the critical energy density e c = E c /N for destabilizing the zone boundary mode vanishes as 1/N 2 , 
independently of dimension. 

C. Large TV limit using the Nonlinear Schrodinger equation 

The large N limit expressions (!■''> .'ill and (|34|l can be derived also by continuum limit considerations. We will derive 
the general expression for any dimension d. The displacement field can be factorized into a complex envelope part ip 
multiplied by the zone boundary mode pattern in d dimensions. 

u nu ..., ni = ^K-^^ St^ ,0+e.a , (35) 

where 



w*,...,* = Vd [4(1 + 3|V| 2 )]. (36) 

Substituting Eq. (|35|l into the FPU lattice equations in d dimensions, a standard procedure [iH leads to the 
following d dimensional Nonlinear Schrodinger (NLS) equation: 

^ + |a^-QV#| 2 = 5 (37) 
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where is the d dimensional Laplacian. The parameters P and Q are derived from the nonlinear dispersion relation 

d 



[4sin 2 ^ + 12M 2 sh 

»=i 



4 It 



(38) 



P 



8 2 uj 



(<?1 =TT,...,q d = 7T, |-0| = 0) 



2%/d 



duj 



(qi = 7T, ...,q d = 7r, IV'I = 0) = -3Vd. 



(39) 
(40) 



Assuming that, at the first stage, modulational instability develops along a single direction a; and that the field 
remains constant along all other directions, one gets the one-dimensional NLS equation 



M Pd 2 i> a 



(41) 



Following the results of the inverse scattering approach any initial distribution of amplitude \if)\ and length A 



ampi 

along x, and constant along all other directions, produces a final localized distribution if |18| 



(l^|A) 2 >^ 2 



(42) 



This means that if the initial state is taken with constant amplitude \ip\ = a on the d-dimensional lattice with N d 
oscillators, the modulational instability threshold is 



(43) 



which coincides with the leading order in Eq. I|33|) . 



III. CHAOTIC BREATHERS 



In this Section, we will discuss what happens when the modulational energy threshold is overcome. The first 
thorough study of this problem can be found in Ref. [2l| , many years after the early pioneering work of Zabusky and 
Deem |4|. Already in Ref. l2lL it has been remarked that an energy localization process takes place, which leads to 
the formation of breathers 26] . This process has been further characterized in terms of time-scales to reach energy 
equipartition and of quantitative localization properties in Ref. The localized structure which emerges after 

modulational instability has been here called "chaotic breather" (CB). The connection between CB formation and 
continuum equations has been discussed in Refs. [2^,[3(J, while the relation with the process of relaxation to energy 
equipartition has been further studied in Ref. |2^. We will briefly recall some features of the localization process in 
one dimension and present new results for two dimensions. 

For long time simulations, we use appropriate symplectic integration schemes in order to preserve as far as possible 
the Hamiltonian structure. For the one dimensional FPU, we adopt a 6th-order Yoshida's algorithm ^3] with a time 
step dt — 0.01; this choice allows us to obtain an energy conservation with a relative accuracy AE/E ranging from 
10 -10 to 10~ 12 . For two dimensions, we use instead the 5-th order symplectic Runge-Kutta-Nystrom algorithm of 
Ref. which gives a similar quality of energy conservation. 

We report in Fig. [2^a) a generic evolution of the one dimensional 7r-mode above the modulation instability critical 
amplitude (a > a c ). The grey scale refers to the energy residing on site n, 

E n = 7^ + ^V(u n+1 - u n ) + ^V(u n - u„_i), (44) 

where the FPU-potential is V(x) = \x 2 + \x A . Figs. E^b), Etc) andEtd) are three successive snapshots of the local 
energy E n along the chain. At short time, a slight modulation of the energy in the system appears (see Fig.Etb)) and 
the 7r-mode is destabilized [l^. Later on, as shown in Fig.Et a )i only a few localized energy packets emerge: they are 
breathers 26] . Inelastic collisions of breathers have a systematic tendency to favour the growth of the big breathers 
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FIG. 3: Time evolution of the local energy l|44jl . In panel (a), the horizontal axis indicates lattice sites and the vertical axis 
is time. The grey scale goes from E n — (white) to the maximum i5 n -value (black). The lower rectangle corresponds to 
< t < 3 10 4 and the upper one to 0.994 10 6 < t < 10 6 . Figs, (b), (c) and (d) show the instantaneous local energy E n along 
the N = 128 chain at three different times. Remark the difference in vertical amplitude in panel (c), when the chaotic breather 
is present. The initial 7r-mode amplitude is a = 0.15 > a c ~ 0.01. 



at the expense of small ones |49|,|50|. Hence, in the course of time, the breather number decreases and only one, of 
very large amplitude, survives (see Fig. Etc)): this is the localized excitation we have called chaotic breather (CB). 
The CB moves along the lattice with an almost ballistic motion: sometimes it stops or reflects. During its motion 
the CB collects energy and its amplitude increases. It is important to note that the CB is never at rest and that it 
propagates with a given subsonic speed 51] . Finally, the CB decays and the system reaches energy equipartition, 
as illustrated in Fig. Eld). The final CB decay is present in all the simulations we have performed, but one cannot 
exclude the existence of examples where the final breather does not disappear. 
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In order to obtain a quantitative characterization of energy localization, we introduce the "participation ratio" 

N 

C (t) = N -if (45) 



which is of order one if Ei — E/N at each site of the chain and of order N if the energy is localized on only one site. 
In Fig. Ufa), Co is reported as a function of time. Initially Co grows, indicating that the energy, evenly distributed 
on the lattice at t = 0, localizes over a few sites. This localized state survives for some time. At later times, Cq 
starts to decrease and finally reaches an asymptotic value Co which is associated with the disappearance of the CB 
(an estimate of Co has been derived in Ref. |27j taking into account energy fluctuations and is reported with a dashed 
line in Fig. Ufa)). At this stage, the energy distribution in Fourier space is flat, i.e. a state of energy equipartition is 
reached. 

As explained in Ref. |2jj , the destruction of the breathers is related to its interaction with low frequency modes 
which arc dominant in the chain after the initial stage. A full study of the scattering of plane waves by FPU-breathers 
would be necessary to quantify this explanation either as anticipated in Ref. |27| or, even better, as recently performed 
by Flach and collaborators for the discrete NLS equation, Klein-Gordon or FPU chains |53l l54j . 

In Fig.^b), we show the finite time largest Lyapunov exponent \i(t) for the same orbit as in Fig.QJa). We observe 
a growth of Ai(t) when the CB emerges on the lattice and a decrease when it begins to dissolve. The peak in Xi(t) 
perfectly coincides with the one in Co- Due to this increase of chaos associated with localization, we have called the 
breather chaotic (although chaos increase could be the result of a more complicated process of interaction with the 
background) . 

In Ref. |27j . the time-scale for the relaxation to equipartition has been found to increase as (E/N)~ 2 in the small 
energy limit. This has been confirmed by the followers of this study |2^, [2^, . Such power law scalings are found 
also for the FPU relaxation starting from long wavelengths jH^I: the so-called FPU problem. We have termed the 
relaxation process which starts from short wavelengths the Anti-FPU problem, just because of the similarities in 
the scaling laws. The main feature of the latter problem is that relaxation to equipartition goes through a complex 
process of localized structures formation well described by breathers or, in the low-amplitude limit, by solitons of the 
NonLinear Schrodinger equation. On the contrary, for the original FPU problem, an initial long wavelength excitation 
breaks up into a train of mKdV-solitons. The final relaxation to equipartition is however due to an energy diffusion 
process which has similar features for both the FPU and the Anti-FPU problem |29|. 

A similar evolution of the local energy 

En.vri 7>^n,m ) + -^V(u n:1n+1 - U n>m ) 

+ 7^(«n-l,m - «n,m) + jV{u n ,m-l ~ U n ,m) ( 46 ) 

is observed for the two-dimensional case (see Fig. |5j). In this figure, we just show the initial evolution which leads 
to breathers formation. As for the one-dimensional case, bigger breathers eat up smaller ones, and finally only one 
breather survives. However, note that, depending on the initial conditions, the simulations do not always lead to the 
coalescence into a single breather, because collisions are more rare in two dimensions than in one. After the formation 
of a few localized structures, one also observes the final relaxation to equipartition, which is not shown in Fig. [SJ This 
latter is instead evident from the time evolution of Co(t), the localization parameter, shown in Fig. EJc): its behavior 
is very similar to the one-dimensional case. Indeed, also the largest finite time Lyapunov exponent behaves similarly 
(see Fig. Eld)). 



IV. SPONTANEOUS LOCALIZATION BY EDGE COOLING 



Breathers play an important role in the non-equilibrium dynamics of the FPU model. The relaxation to equiparti- 
tion of the zone-boundary mode, analyzed thus far in this paper, is one example. Another interesting case in which 
breathers emerge spontaneously is when the lattice is cooled at its boundaries [^J l32l l33l l34l l35| . This process may 
be thought of as modeling a non-equilibrium process where energy exchange in the bulk is much slower than at the 
surface. Although in this case the dynamics is dissipative, it turns out that there is a deep connection with the original 
Hamiltonian model. 
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FIG. 4: Panel (a) presents the evolution of Co(t) of formula l|45[l for the one-dimensional FPU lattice with N = 128 oscillators, 
initialized on the 7r-mode with an amplitude a — 0.126 > a c ~ 0.010. The dashed line indicates the equilibrium value Co = 1.795. 
Panel (b) presents the corresponding finite time largest Lyapunov exponent. Panel (c) shows Co(t) for the two-dimensional 
FPU lattice with 20*20 oscillators, initialized on the (-7T, 7r)-mode with an amplitude a — 0.425 > a c ~ 0.045. Panel (d) presents 
the finite time largest Lyapunov exponent for two dimensions. 



Specifically, when modelling this process, we add a dissipative term —r]u n in the r.h.s. of Eq. (JTJ when n = 1 and 
n = N. Similarly, in two dimensions, the same term is introduced at all edge sites. The parameter rj controls the 
strength of the coupling with the external reservoirs at zero temperature. Since we are interested in the exchange 
of energy between a finite system and the environment, we shall consider either free-ends or fixed-ends boundary 
conditions. 

In a typical simulation, first an equilibrium micro-state is generated by letting the Hamiltonian system (r; = 0) 
evolve for a sufficiently long transient. Then, the dissipative dynamics (77 > 0) is started. The initial condition for the 
Hamiltonian transient is assigned by setting all relative displacements to zero and by drawing velocities at random 
from a Gaussian distribution. The velocities are then rescaled by a suitable factor to fix the desired value of the initial 
energy E(0) (see Ref. [33| for the details). 

The one-dimensional numerical simulations reveal that the dissipation rate of the energy is dominated by two 
sequential effects, that characterize the pathway to localization. In the first stage of the energy release process, 
the relaxation law undergoes a crossover from the exponential exp[— t/r ] to the power law (t/r ) _1 / 2 , where t = 
N/(2rf) sets the shortest time scale of the system (see Fig. |SJ) |34). Asymptotically, energy reaches a plateau and, 
correspondingly, the localization parameter Co(i) also saturates (see the inset of Fig. |SJ). 

The crossover at t is a signature of the hierarchical nature of the early stages of the process. In the harmonic 
approximation, if one adds a small dissipation, the frequency uj of each linear mode acquires an imaginary part y(u>), 
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FIG. 5: Local energy l|46|l surface plots for the two-dimensional FPU lattice with 20*20 oscillators, initialized on the (it, 7r)-mode 
with an amplitude a — 0.225 > a c ~ 0.045. Snapshots at four different times t are shown. Breathers form after a coalescence 
process similarly to the one-dimensional case. The mobility of the breathers is evident and one also observes in the last panel 
the final decrease. 

which represents its damping rate. Initially, it is only the fastest mode which determines the energy relaxation rate. 
As time goes on, past t ~ r , it is the full spectrum of decay times of the linear modes that sets the rules of energy 
relaxation. Actually, it turns out that at this stage, for small nonlinearities, the system behaves approximately as its 
linear counterpart (see Fig.|B). In particular, a perturbative calculation to first order in 7 confirms that |34| 



where the density of states g{u>) is derived from the dispersion relation u)(q) = 2 sin(g/2), with q = -nk/N for free-ends 
boundary conditions, and q = ir(k + 1)/(N + 1) for fixed-ends boundary conditions (k = 0, . . . N — 1). 

The time range after the crossover coincides with the onset of localization. Now the dynamics is significantly 
affected by the spontaneous appearance of breathers. As the latter exhibit a very weak interaction among themselves 
and with the boundaries, the energy release undergoes a slowing down, thus freezing the system in a quasi-stationary 
configuration far from thermal equilibrium. 

Such residual state is characterized by the presence of one, highly energetic localized object (possibly accompanied 
by a few much smaller ones), which is mobile and alternates periods of rest and erratic motion, as shown in Fig. 
This behaviour is reminiscent of the chaotic breathers which emerge in the Hamiltonian system from modulational 
instability of band-edge modes, discussed in Section ITTTI fsee also Fig. |3J. 

But what do we know of the mechanisms leading to spontaneous localization in the presence of dissipation? As 
noticed above, there is evidence that it is the modulation instability of short lattice waves that triggers the formation 
of localized structures. This hypothesis is supported by the observation that the emergence of spatial patterns in the 
early stages of the relaxation is intimately related to how dissipation acts on vibrational modes of different wavelength. 
In particular, if the system is swiftly enough depleted of long-wavelength modes, the instability associated with the 
band-edge waves may effectively trigger the process of localization. A key test for the above hypothesis is offered by 
the nature of the boundary conditions. In the case of free-ends, the modes of small wavenumber indeed disappear very 
fast, in conjunction with the onset of spontaneous localization. On the contrary, in the presence of fixed boundaries 
the former turn out to be as long-lived as the modes with large wavenumbers. The corresponding numerical evidence 
is that hardly no localization is observed in this case. Rather, the energy decays following exactly the behavior of the 




for t < t c 
for t > T c 



(47) 
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FIG. 6: Log-log plot of — log[E(t)/E(0)] versus time for an FPU chain with free-ends boundary conditions. In this representa- 
tion, an exponential is a straight line with slope one. Symbols arc the results of numerical simulations averaged over 20 initial 
conditions. The dashed line is a plot of the theoretical result I47H 33] for a harmonic chain. The arrow indicates the crossover 
time to. Parameters are N — 100 and r\ = 0.1. 
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FIG. 7: Space-time contour plot of the site energies l)44|l . Time flows up and the horizontal axis is the site index. Parameters 
are N = 100, rj = 0.1 and E(0)/N = 1. 



harmonic chain l|47|) 33]. This scenario can be observed directly, by computing the time-dependent spatial power 
spectrum S(q,t) during the relaxation (see Fig. 0. 

It is possible to get a quantitative confirmation of the above hypothesis by calculating the exact relaxation spectrum 
7(w) of the linearized system. In the harmonic approximation, the equations of motion may be written in matrix 
form as 



U = AU, 



(48) 
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FIG. 8: Surface plot of the time-dependent spatial spectrum of particle velocities for the one-dimensional FPU lattice, (a) 
Fixed-ends boundary conditions, (b) Free-ends boundary conditions. Parameters are N = 100, r\ = 0.1 and E(0)/N = 1. 



where U = (ui, 



, un) t is a 2N column vector, and A is the matrix 
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(49) 



The tri-diagonal matrix of force constants K nm also contains the information on the type of boundary conditions, 
whereas the matrix B nm — S nm (S m i + <5 m 7v) describes the coupling with the environment. The spectrum of damping 
rates can be calculated by straight diagonalization of matrix A. In Fig. [5] we plot 7(0;), where 7 is the opposite of the 
real part of an eigenvalue of A and lu its imaginary part. This representation of the relaxation spectrum is preferable 
with respect to the one in terms of wave-numbers, since the vibration frequencies are shifted as a consequence of the 
damping. 

This calculation confirms that the free-ends and fixed-ends systems display considerably different behaviors. In the 
former case, the least damped modes are the short-wavelength ones (cu « 2, the band-edge frequency), the smallest 
damping constant being 7(2) w ir 2 rj/2N 3 , while the most damped modes are the ones in the vicinity of u> ~ 1/N. On 
the contrary, for fixed ends, the most damped modes are those around the band center (uj ~ y/2) while both short- 
and long- wavelength waves dissipate very weakly (7 <~ 1/N 3 ). 

The above analysis helps understanding why spontaneous localization is strongly inhibited if the system is trapped 
between rigid walls, thus in parallel unveiling the role of modulational instability in the process. 

We have also performed similar numerical simulations for the two-dimensional FPU model (|19J) . with dissipation 
added at the edges and free-ends boundary conditions. Remarkably, the asymptotic scenario changes. The quasi- 
stationary state is now a static collection of tightly-packed localized objects, arranged in a sort of random lattice (see 
Fig. 110(1 . Moreover, it turns out that spontaneous localization in two dimensions is a thermally-activated phenomenon, 
described by an Arrhenius law for the average breather density, where the parameter that controls the strength of 
thermal fluctuations is the initial energy density E(0)/N (34|. The origin of this behaviour is that, in the two- 
dimensional FPU system, discrete breathers may be excited only above a certain energy threshold (2^, as discussed 
in Section Til Bl Despite the different nature of the asymptotic state, the onset of localization follows the same path 
in one and two dimensions [34|. well described the hierarchy of relaxation times underlying Eq. (|47H . In particular, a 
crossover is observed from the exponential exp[— 2t/r ] to the power law (i/r ) _1 . 



V. CONCLUSIONS 



In this paper, we have presented a detailed analysis of the zone-boundary mode modulational instability for the 
FPU lattice in both one and higher dimensions. Formulas for the critical amplitude have been derived analytically 
and compare very well with numerics for all system sizes. We have extended to two dimensions the study of the 
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FIG. 9: Damping rates 7(o>) for an harmonic chain with N = 100 and r\ — 0.1. Free-ends boundary conditions (stars) and 
fixed-ends boundary conditions (diamonds). 
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FIG. 10: 2D FPU lattice, site energies in the residual state state. Parameters are: N — 80, r\ = 0.1, E(0)/N — 1. 

process which leads to the formation of chaotic breathers. The physical picture is similar to the one-dimensional case. 
Besides that, we make the bridge between breathers created by modulational instability of plane waves and the ones 
formed when extracting energy from the boundaries: similarities and differences are highlighted. 

All results on modulational instability of zone-boundary modes can be straightforwardly extended to other initial 
modes and, correspondingly, instability rates can be derived. This has already partially been done in Ref. |14| and 
compares very well with the numerical results by Yoshimura |55| . This author has recently reanalyzed the problem j5|J 
to determine the growth rates for generic nonlincaritics in the high energy region, obtaining exact results based on 
Mathicu's equation. 

Going to many-modes initial excitations, it has been remarked that instability thresholds depend on relative 
amplitudes and not only on the total energy . Although this makes the study of the problem extremely involved, 
we believe that a detailed study of some selected group of modes, which play some special role in FPU dynamics, 
could be interesting. The method discussed in this paper could be adapted to treat this problem. Historically, the 
first study is in the paper by Bivins, Metropolis and Pasta himself [53, where the authors tackle the problem by 
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studying numerically the instabilities of coupled Mathieu's equations. 

The study we have presented in this paper of the two-dimensional FPU lattice is extremely preliminary and further 
analyses are needed. In particular, the full process of relaxation to energy equipartion and the associated time scales 
have not been studied at all. Preliminary results on the relaxation process in two dimensions from low frequency 
initial states seem to indicate a faster evolution to equipartition |58j . A similar analysis for high frequencies remains 
to be performed. 

In one-dimensional stud ies, a connection between the average modulation instability rates and the Lyapunov expo- 
nents has been suggested [ill EjL, Recently [5!j , high frequency exact solutions have been used in the context of a 
differential geometric approach [6(j to obtain accurate estimates of the largest Lyapunov exponent. Similar studies 
could be performed for the two-dimensional FPU lattice and the corresponding scaling laws with respect to energy 
density could be obtained. 

The study we have reported in the last Section about lattices that are cooled at the boundaries points out the 
similarity of the localized objects obtained in the long time limit with chaotic breathers. However, this resemblance, 
although convincing, is only qualitative. Quantitative studies on the comparison of these breathers with the chaotic 
ones obtained from modulational instability should be performed. 
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